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We compute the spectra and fluxes of the velocity and temperature fields in Rayleigh-Benard 
convection in turbulent regime for a wide range of Prandtl numbers using pseudo-spectral simula- 



tions on 512 3 grids. Our spectral and flux results support the Kolmogorov-Obukhov (KO) scaling 
for zero Prandtl number and low Prandtl number (P = 0.02) convection. The KO scaling for 
the velocity field in zero-Prandtl number and low-Prandtl number convection is because of the 
weak buoyancy in the inertial range (buoyancy is active only at the very low wavenumbers). We 
also observe that for intermediate Prandtl numbers (P = 0.2) the KO scaling fits better with the 
numerical results than the Bolgiano-Obukhov (BO) scaling. For large Prandtl number (P = 6.8), 
the spectra and flux results are somewhat inconclusive on the validity of the KO or BO scaling, yet 

m 

| the BO scaling is preferred over the KO scaling for these cases. The numerical results for P = 1 is 

rather inconclusive. 
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I. INTRODUCTION 

Turbulent convection is one of the most challenging problems of classical physics A 
large number of work on convection have been done for an idealized version called Rayleigh- 
Benard convection (RBC) in which the fluid is heated between two parallel plates. The 
convective flow properties depend on two nondimensional parameters: the Rayleigh number 
(proportional to the buoyancy force) and the Prandtl number (the ratio of kinematic viscosity 
and thermal diffusivity) . The convective flow becomes turbulent when the Rayleigh number 
is much larger than the critical Rayleigh number. One of the important topics in the study 
of convective turbulence is the scaling of energy spectra and energy fluxes of the velocity 
and temperature fields in the inertial range. In this paper we compute these quantities using 
direct numerical simulation (DNS) and compare them with the predicted values from the 
existing phenomenologies. 

The energy spectra and fluxes for convective turbulence are more complex than those for 

nn 

fluid turbulence due to the presence of the buoyancy force [2|, 13J. For stable stratified fluid 
convection, Bolgiano 4] and Obukhov jl5] proposed dual cascade in the inertial range. For 
small wavenumbers (large length scale), they predicted dominance of the buoyancy force 
over the inertial force leading to the velocity and temperature spectra as k~ lx ^ and k~ 7 ^ 5 
respectively, where k is the wavenumber. In this regime, the energy flux of the temperature 
field is constant, while the flux of the velocity field varies as k~ 4 ^ 5 . For the intermediate 

n n 

wavenumbers, Bolgiano [4J and Obukhov [5[ conjectured dominance of the inertial force 
over the buoyancy force. Consequently the temperature field evolves as a passive scalar, and 
both the velocity and temperature fields have Kolmogorov's energy spectrum (A; -5 / 3 ) and 
constant energy fluxes ^-6|. The length scale that separates these two different regimes of 
energy cascades is called the "Bolgiano length" (/#). 



Later Procaccia and Zeitak |7j], L'vov 



and Falkovich and L'vov {9] proposed the 



same scaling for Rayleigh- Benard convection. In convective turbulence, for scales above the 
Bolgiano length (I > /#), the kinetic energy spectrum (E u (k)) and the entropy spectrum 
(E e (k)) follow the Bolgiano- Obukhov (BO) scaling 



E u (k) = C k (e°)Hag)r>k--, 
E\k) = C e (e 9 )l(agrik- 7 s, 
W{k) = C f {e )Hag)h-l 



(1) 
(2) 
(3) 



and for I < Ib, spectra follow Kolmogorov-Obukhov (KO) scaling 



E u {k) = K ko {e u )*k-^ 
E\k) = K e e e (e u )-h-l 



(4) 
(5) 



where IT" is the kinetic energy flux, e u and e e are the kinetic and entropy dissipation rates 
respectively, a is the thermal expansion coefficient of the fluid, and g is the acceleration due 
to gravity. Note that in literature, the spectrum and the flux of the temperature field are 
also referred to as the "entropy spectrum" and "entropy flux" respectively. 

The Bolgiano length Is has the following dependence on the convective parameters: 
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Nu?d 



(6) 



where Nu is the Nusselt number (dimensionless heat flux), R is the Rayleigh number, P is 
the Prandtl number, and d is the vertical height of the container. Grossmann and L'vov [lo| 
and Cioni et al. HI] argued that for P < 1, Bolgiano length is of the order of container's size. 
Hence, only KO scaling is expected in the inertial regime for low Prandtl number (low-P) 
convection. For large-Prandtl number (large-P) convection, Ib lies in the inertial regime, 
hence mixed scaling is expected. Several exact relationships connecting e u , e 6 , Nu, R, and P 
have been derived for homogeneous convective turbulence. Shraiman and Siggia 12J derived 
that 
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Researchers have attempted to test the above scaling predictions [Eqs. (1-5)] using ex- 
periments and numerical simulations (to be described later in this section). Yet, the scaling 
of convective turbulence has not been conclusively established. In a recent review, Lohse 
and Xia s| described these results critically and exhaustively. The inconsistencies of the 
scaling predictions with numerical and experimental results are attributed to the drastic 
assumptions made in the scaling arguments. In the theory described above both thermal 
and viscous boundary layers are not considered appropriately. Shraiman and Siggia 12] and 
Grossmann and Lohse 13[ showed that the properties of the fluctuations in the boundary 
layer and in the bulk are rather different. Experiments and numerical simulations reveal 



that the fields in the boundary layer are highly inhomogeneous and anisotropic, while the 
bulk flow is somewhat homogeneous and isotropic. Hence the above scaling arguments are 
expected to hold only in the bulk, if at all. The computation of the Bolgiano length Ib 
[Eq. (|6l) | assumes uniform dissipation, which is not valid in the boundary layer. Calzavarini 
et al. [14| have computed Ib for different layers in the convective fluid; they report that ls/d 
is small near the walls (in the boundary layer), but ls/d ~ 1 in the bulk. In brief, the pres- 
ence of boundary layers, a single Is for the whole fluid, inhomogeneity and anisotropy of the 
flow are some of the features that possibly make the above scaling arguments inconsistent 
with realistic experiments and simulations 

To disentangle various complexities mentioned above, some researchers have idealized the 
geometry of RBC even further. For example Borue and Orszag [l5], Skandera et al. 16] 
considered convection in a periodic box (with thermal gradients along the vertical) and 
obtained KO scaling. This feature removes the effects of the viscous and thermal boundary 



layers on the bu 
boundary layer 



k, and hints that BO scaling is possibly due to the thermal forcing in the 
3j. In the present paper we consider free-slip and conducting boundary 
conditions in which viscous boundary layer is insignificant, while the thermal boundary 
layer is present. We expect that our numerical results will suppress the effects of viscous 
boundary layers, and may possibly provide scaling for the bulk convective flow. 

In the following discussion we briefly review the experimental studies that attempt to 
test the above phenomenology of RBC. Many convection experiments measured the veloc- 
ity and temperature fields only at fixed locations of the apparatus. For such experiments 
"Local Taylor hypothesis" is invoked to relate the frequency spectrum to the wavenumber 



spectrum 



171 ]. However, in some experiments, high resolution spatial velocity and tem- 



perature fields have been measured for computing the above mentioned spectra; experiments 
by Mashiko et al. 18[ and Sun et al. 19[ belong to this category of experiments. Chilla et 
al. 2(3], Zhou and Xia 21], and Shang and Xia 22] carried out convection experiments on 
water (P ps 7) at large Rayleigh number and found the energy spectrum to be consistent 



with BO scaling. Heslot et al. [23| and Castaing 24| measured frequency power spectrum of 
the temperature field in He gas (0.65 < P < 1.5) and found the spectrum to be consistent 



with KO scaling. Wu et al. 25] however reported BO scaling for Helium gas through fre- 



quenc y sp ectrum measurements of temperature. Ashkenazi and Steinberg 26j] and Mashiko 
et al. 18| performed convection experiments for SF 6 (1 < P < 93) and mercury respectively 



and reported the BO scaling for them. Niemela et al. [27j measured temperature time series 



in He gas and reported presence of both the KO and BO scaling. Cioni et al. S 



carried out 



experiments on mercury (P ~ 0.02, a low-P fluid), and reported KO frequency spectrum 
for it. Thus the outcome of these experiments are somewhat inconclusive on the validity of 
the phenomenologies for RBC, yet majority appear to support the BO scaling for large-P 
convection, and the KO scaling for low-P scaling. 

Numerical experiments provide important clues in the study of turbulence. A series of nu- 
merical simulations of RBC have been performed to test the KO and BO scaling. Grossmann 



and Lohse kill simulated RB fluid with P = 1 under Fourier- Weierstrass approximation 
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and reported KO scaling. Borue and Orszag [15[ and Skandera et al. [16j performed pseudo- 
spectral simulation on P = 1 fluid with periodic boundary conditions on all directions and 
found consistency with KO scaling. Vincent and Yuen 31| performed spectral simulation for 
P = 1 and R = 10 8 using free-slip boundary conditions and reported —5/3 and —3 spectral 
indices for the temperature and velocity fields respectively. They however find dual branches 



in the entropy spectrum. Paul et al. 



32| also observed dual entropy spectrum in their 2D 



spectral simulations with free-slip boundary conditions, albeit at lower Rayleigh numbers. 
Rincon 33J performed a numerical simulation using higher order finite-difference scheme to 



study the effects of inhomogeneity and anisotropy on the scaling of the energy spectra; for 
Ra = 10 6 and P = 1 on 256 x 256 x 128 grids with free-slip boundary conditions, he reported 
that the numerical results are inconclusive in identifying a definite spectral slope. Kerr 34] 
used pseudo-spectral method for his simulations of P = 0.7 fluid ( air) under no-slip bound- 
ary conditions and observed KO scaling. Camussi and Verzicco 35] performed numerical 
simulations for cylindrical geometry using finite difference method; they found both velocity 
and temperature spectral exponents to be —7/5, which is inconsistent with both the KO 
and BO scaling. They attribute this anomaly to inhomogeneities and anisotropy of the flow 
near the boundaries. On the whole, numerical results indicate uncertainty in the tests of 
the convective phenomenology. 

Another way to investigate turbulent scaling is through the structure function calcula- 
tions. Following Kolmogorov, Yakhot 36[ derived an exact analytical form for the third 
order structure function for the BO scaling. Sun et al. 19[ computed the structure function 
of the velocity and the temperature fields using the data of their convection experiments 
on water and reported the KO scaling. Kunnen et al. 3t| performed similar calculations 



for Helium gas (both experiments and numerical simulation) and observed the BO scaling. 
Calzavarini et al. [l^] computed third-order structure function using Lattice Boltzmann 
method for P = 1 and reported the BO scaling. Hence structure function studies too are 
inconclusive on the validity of the BO or KO scaling. 

In this paper we compute the energy spectra and cascade rates for the velocity and 
temperature fields using pseudo-spectral method on 512 3 grids with free-slip boundary con- 
ditions. Our computations include zero-Prandtl number (zero-P), low-P, and large-P con- 
vection regimes (P = 0,0.02,0.2,1,6.8), hence we have reasonable number of numerical 
runs to test the convective turbulence phenomenology. We also construct phenomenological 
arguments to understand zero-P and low-P numerical results. 

The outline of the paper is as follows. Section II contains the dynamical equations and 
the definitions of the energy spectra and fluxes. The details and results of our numerical 
simulations are discussed in section III. We conclude in Sec. IV. 



II. GOVERNING EQUATIONS 

We numerically solve the nondimensionalized Rayleigh-Benard equations under the 



Boussinesq approximation 
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^ + (u ■ V)u = -Va + R9z + J ^V 2 u, (9) 

P + (u ■ V)*) =u 3 + y|v 2 fl, (10) 
Vu = (11) 

where u = (ui,u 2 ,u 3 ) is the velocity field, 9 is the perturbations in the temperature field 
from the mean temperature, a is the deviation of pressure from the conduction state, R = 
ag(AT)d 3 / vk is the Rayleigh number, P = v/nis the Prandtl number, and z is the buoyancy 
direction. Here v and k are the kinematic viscosity and thermal diffusivity respectively, d is 
the vertical height of the container, and AT is the temperature difference between the plates. 
For the nondimensionalization we have used d as the length scale, ^a(AT)gd as the velocity 
scale, and z/( AT) /k as the temperature scale. For large-P convection, the temperature scale 
is taken as AT, and the governing equations are altered accordingly. 

Zero-Prandtl number (Zero-P) convection is the limiting case of low-P convection. The 



corresponding dimensionless equations for zero-P convection are 

^ + (u • V)u = -Vo- + BB + V 2 u, (12) 

u 3 + v 2 e = o (13) 

Here we use d as the length scale, v/d as the velocity scale, and v(AT)/k as the temperature 
scale. 

Boundary conditions of the systems strongly affect the properties of the convective 



flow 



13 



301 ] . We employ free-slip and conducting boundary conditions on the horizontal 



plates, hence 

^3 = d 3 ux = d 3 U2 = 8 = 0, at z = 0, 1. (14) 

Periodic boundary conditions are applied along the horizontal directions. Consequently the 
velocity and temperature fields are expanded in terms of basis functions as 



y*][ui,2{ix, i y , 0) + ^2 ui j2 (i x , i y , i z )2 cos(k z z)] exp i(k x x + k y y) (15) 



u 3 (x,y,z) = u 3 (i x , i y , i z )2 sin(k z z) exp i(k x x + k y y) (16) 
9(x,y,z) = y j 6(i x , i y , i z )2 sm(k z z) exp i(k x x + k y y) (17) 



where (i x ,i y ,iz) are the grid indices with k x = i x ir/\/2, k y = iyir/y/2, and k 



nir. 



The energy spectra of the velocity field (E u (k)) and the temperature fields (E e (k)) are 
defined as 



E u (k)= Yl \n k ')\ 2 > ( 18 ) 

k<k'<k+l 

E\k)= ^( k ')| 2 - (19) 

k<k'<k+l 

Here the sum is being performed over the Fourier modes in the shell [k,k + 1). We will 
compute these spectra numerically at the steady state. Note that the magnitude of the 
wavevectors in Fourier space is 

1 /2 

k = (M/v^) 2 + M/v^f + M) 2 ) • (20) 
Here we use the fact that the aspect ratio of the box is 2a/2 [l|. 



The energy flux is a measure of the nonlinear energy transfers in turbulence 39l-l41|. 
The energy flux for a given wavenumber sphere is the total energy transferred from the 
modes within the sphere to the modes outside the sphere. The energy flux for fluid and 
magnetohydrodynamic turbulence has been studied in great detail. However there are only 
a small number of work on the flux computations in convective turbulence 15|, |l6|, |42(. Toh 



and Suzuki |42| defined the kinetic energy flux U u (ko) and the entropy flux n e (/c ) based on 
Kraichnan formalism 41] as 

(21) 



n u (fc 



k>ko p,q<ko 

m (p)it n (q) 

n e (M = \ J2 E 4w(k-u( q )) 

k>ko p,q<ko 

x(r(k)^( P )) 



(22) 



These quantities represent the net cascade of \0\ 2 /2 and |u| 2 /2 respectively from the modes 
within the wavenumber sphere of radius ko to the modes outside of the sphere. 

The energy fluxes defined above can also be defined quite conveniently using the "mode- 
to-mode energy transfers" formalism discussed in Verma Q|. According to this formalism, 
the kinetic energy flux and the entropy flux are 



n ^o) = ^^4, P+q ^([k-u(q)][ U *(k)-u(p)]) 

k>ko p<ko 

n e (M = EEW(P^«(q)]W(p)]) 



(23) 
(24) 



k>ko p<k(, 

where 9" represents the imaginary part of the argument. We compute the spectra and fluxes 
of the velocity and temperature fields using numerical simulations {40]. These results will 
be described in the next section. 



III. NUMERICAL SIMULATIONS AND RESULTS 

As described in the previous section, the dynamical equations of RBC are Eqs. (IMTTT) 
for low-P convection and Eqs. fl 1 1 til 3 f) for zero-P convection. The equations for large-P con- 
vection are similar. We solve these equations numerically using a pseudo-spectral method 
under free-slip boundary conditions for the horizontal plates, and periodic boundary condi- 



tions along the horizontal directions. The expansion of the velocity and temperature fields 
are given in Eqs. ( TTBTTPTj) . 

The unidirectional initial energy and entropy spectra for the initial conditions are of the 
form: 

m = (k4 k qr+" eM ~ bk1 ' 11 (25) 



where b = 0.02, q = 1.5, a = 2.8/12, and a as a free parameter 43j. The initial phases are 
generated randomly. Time stepping of dynamical equations are carried out using fourth- 
order Runge-Kutta (RK4) scheme. We start our simulation on a smaller grid and run it 
until the steady state is reached. We then use the steady solution of the lower grid as an 
initial condition for simulations on larger grid size at a larger R. We continue this procedure 
until turbulence state is reached. The final runs were performed on 512 3 grid for 20 large 
eddy turnover time on 8 nodes and 16 nodes of EKA, the supercomputer at Computational 
Research Laboratory, Pune. Zero-P convection runs were performed on 256 3 grid. The 
kmax^i where rj is the Kolmogorov length, for our simulations are always greater than one 
indicating that our simulations are well resolved. 

For the energy flux calculation, we divide the wavenumber space into 20 shells. The first 
three shells are k = (0,2), [2,4), and [4,8), and the last shell contains all modes beyond 
k = 568. Between k = 8 and k = 568, the wavenumber space is split into shells bounded by 
[k n , k n+l ) with k n = 8x 2< n ~^ where s = (1/15) ln 2 (568/8). 

For free-slip boundary conditions, the viscous boundary layer is practically absent, while 
the thermal boundary layer is significant 44]. To probe the existence of these boundary 
layers we compute the average value of the rms velocity fluctuations and the temperature 
field over the horizontal planes. Figures DJa) and H^b) exhibit these quantities as a function 
of vertical height for R = 6.6 x 10 6 and P = 6.8. We observe a thin thermal boundary 
layer near the horizontal plates (5/d ~ 0.05). The slow variation in the velocity fluctuations 
however demonstrates the insignificance of the viscous boundary layer. Our results are 
consistent with earlier work on boundary layers 441 ] . 

No-slip boundary conditions are encountered more often in convection experiments. In 
our paper we are using free-slip boundary conditions for simplification. An added advantage 
of the free-slip boundary conditions could be a reduction of the complexity of the viscous 
boundary layer; as a result, the energy spectrum of the flow may reflect the bulk properties. 



Thus we may be able to probe the validity of the KO or BO scaling for the bulk flow using 
these simulations. Note that several properties of the convection are the same for both 
free-slip and no-slip boundary conditions, e.g., the scaling exponent of the Nusselt number 



45]. 



vs. Rayleigh number is the same for the two boundary conditions 

We choose five representative Prandtl numbers P = 0,0.02,0.2,1,6.8 for our energy 
spectra and flux studies. We compute energy spectra and fluxes, and Nusselt number using 
the numerically generated data. We also compute e u , e e using the exact relationships [Eqs. (13- 
IHJ)]. Kolmogorov's dissipation wavenumber (kd) and "Kolmogorov's diffusion wavenumber" 
(k c ) are also computed using the phenomenology of passive scalar turbulence 391 ]: 

1/4 



< 26 > 

K = P 3/4 (2g) 
kd 

In Table H] we list the numerically computed and the estimated e u and e e , k c , kd, and inverse of 
the Bolgiano length. The estimated values of e u and e e match quite well with the simulation 
results, thus validating our simulations. 

In the following we will discuss our numerical results on the energy spectra and fluxes for 
various Prandtl numbers. 



A. Prandtl number P = 

For P = 0, the temperature fluctuations can be expressed as 0(k) = M 3 (k)/fc 2 [see 
Eq. (j!3p ]. Consequently E e (k) « E u (k)/k i . Hence the entropy spectrum is very steep 
for zero-P convection, and we can safely assume that the velocity field is buoyantly forced 
only at very large scales (small k). Hence, Kolmogorov's argument for the fluid turbulence 
must be valid for zero-Prandtl number convection. These arguments closely resemble the 



mathematical derivation of Spiegel 46]. 



We performed DNS for P = at R = 1.97 x 10 4 and computed the energy spectrum using 
the steady-state data. In Fig. [2] we plot the compensated energy spectra E(k)k 5 ^ 3 (KO) and 
E(k)k u ^ 5 (BO). Clearly the numerical plots fit better with the KO scaling than the BO 
scaling, thus verifying the above phenomenological arguments. Using the simulation data 



we also compute the kinetic energy flux that is plotted in Fig. |3j The kinetic energy flux is 
flat in the inertial range, in agreement with the KO scaling. The Kolmogorov constant for 
P = is around 1.8 (with the significant errors) which is in a reasonable agreement with 
the expected value of 1.6 (Kolmogorov's constant for the fluid turbulence). 

In the next subsection we will discuss the numerical results for P = 0.02 that can be a 
representative case for low-Prandtl number convection. 



B. Prandtl number P = 0.02 



In the previous subsection we showed that Kolmogorov's scaling (KO) is expected to hold 
for zero-P convection because buoyancy for this case is dominant at very small wavenumbers. 
Here we will attempt to extend the above arguments to low-Prandtl number convection. The 
inertial range for the velocity and temperature fields extends almost up to the Kolmogorov 
dissipative wavenumber (kd) and the Kolmogorov diffusive wavenumber (k c ) respectively. For 
low-P convection, where thermal diffusivity dominates kinematic viscosity, we expect k c <C 
kd [see Eq. ([281]. According to the turbulence phenomenology of passive scalar turbulence, 
E (k) is a power law for k < k c , and it decays exponentially for k > k c . Hence the buoyancy, 
which is proportional to 9 (cf. Eq. ([9])), would be dominant only for low wavenumbers 
(k < k c ), and we expect Kolmogorov's spectrum for the kinetic energy E u (k) for k c < k < k^. 

According to Eq. (f2~7j) . for small e u and large k, k c could be rather small. Under such 
situations the above phenomenological arguments indicate that the velocity field follows Kol- 
mogorov's spectrum, and the temperature field has diffusive energy spectrum. Interestingly, 
the above arguments for low-Prandtl number convections are consistent with the zen> p 



convection for which k c — > (asymptotic case). As argued by Grossmann and L'vov 10], 



the Bolgiano length for low-P convection could be of the order of the box size, so the BO 
scaling is not expected for low-P convection. 

When k c is large, we need more rigorous theoretical arguments to predict the energy 
spectra for k < k c . Possibly, the buoyancy term is irrelevant in "renormalization group" 
sense (see 0|), and both the velocity and temperature fields may follow the KO scaling for 
k < k c in the inertial range. This scenario is observed for P = 0.2 that will be discussed in 
the next subsection. 

In the following discussions we will compare the above phenomenological predictions with 



numerical results. For P = 0.02 we perform numerical simulation at R = 2.6 x 10 6 which is 
at the lower end of turbulent convection regime. We compute k c , ka, e u , e e , lg , and energy 
spectra and fluxes using the numerical data. As evident from the entries of Table [J k c ~ 25 
which is much smaller than kd — 470. According to the arguments given above, we expect 
a diffusive entropy spectrum for k c < k < kd- We do not expect to observe the KO scaling 
for k < k c since k c is too small. 

Figure H] contains the compensated kinetic energy spectra for the KO and BO scaling. The 
KO scaling fits better with the numerical data than the BO scaling, consistent with the above 
phenomenological arguments. Figure |5] exhibits entropy spectrum that contains two distinct 
branches similar to that observed by Vincent and Yuen 3l| and Paul et al. 32[ in their 2D 
spectral simulations with similar boundary conditions as ours. In Appendix A we construct 
phenomenological arguments based on energy equations and numerical results to estimate 
the values of the temperature modes 9(0, 0,2n). We observe that the maximum entropy 
transfers from the modes 9(n, 0,n) and 6(0, n,n) are to the modes #(0,0, 2n) (the three 
indices are i x , i y , and % z respectively). These arguments lead to predictions that 8(0, 0, 2n) ~ 
-1/(2twt) and E e (2n) ~ l/(4n 2 7r 2 ). For P = 0.02 we have listed the values of 6(0, 0, 2n) for 
n = 1 to 4 in Table [III Here #(0,0,2) ~ -1/2tt, but for higher n's, |0(O,O,2ra)| < l/(2wr), 
possibly due to significant entropy transfers to other modes or due to higher thermal diffusion 
for low-P convection. As we will show later, the relationshop 9(0,0, 2n) — — l/(2mr) works 
quite well for large-P convection. However, a common feature borne out for all Prandtl 
number is that the entropy contents of 6(0, 0, 2n) modes are much larger that the other 
thermal Fourier modes, consequently yielding two branches of entropy spectrum. 

We compare the entropy spectrum with both power law and exponential fits (see Fig- 
ure[5]). As evident from the figure, E e (k) ~ exp(— ak) (the inset), which is in agreement with 
the phenomenological arguments given in the beginning of the subsection. We complement 
our spectral analysis with energy flux studies. Figure [6] shows the kinetic energy and en- 
tropy fluxes. The kinetic energy flux is flat for more than a decade indicating Kolmogorov's 
spectrum for the velocity field, in agreement with the KO scaling for the velocity field. The 
entropy flux however drops sharply, consistent with the exponential nature of the entropy 
spectrum. Using Eqs. ( ITllHi) we compute e u and e e that are quite close to the numerically 
computed energy and entropy fluxes (see Tabled]). Also, the numerical estimate of k c and 
kd are in general agreement with the spectra and flux plots. 



On the whole, the numerical results for P = 0.02, which is a representative of low-P 
convection, appear to favour KO scaling for the velocity field. The temperature spectrum 
appears to be diffusive for the most wavenumber region. These numerical results are in good 
agreement with the phenomenological arguments presented above for low-P convection. 

In the next subsection we report energy spectra and fluxes for P = 0.2. 

C. Prandtl number P = 0.2 

Next we present our numerical results for P = 0.2 at R = 6.6 x 10 6 . In Fig. [7] we plot 
the compensated kinetic energy spectra E u (k)k 5/3 (KO) and E u (k)k u ^ 5 (BO). Even though 
both the BO and the KO scaling do not fit very well with the numerically computed energy 
spectrum, yet the KO scaling is in better agreement with the numerical data than the BO 
scaling. 

In Fig. [8] we plot the entropy spectrum, which has significant inertial range. Note that 
k c = 68 (see Table H]). We obtain bi-spectra similar to that for P = 0.02. As described in 
the Appendix A, the upper curves represent the spectrum of the Fourier modes #(0,0, 2n), 
and it matches reasonably well with k~ 2 spectrum. As evident from the entries of Table HIl 
8(0, 0, 2n) matches with —l/(2mr) within a factor of 2. The lower curve however appears to 
fit better with the KO scaling than the BO scaling. Note that the upper curve of the entropy 
spectrum contains small number of Fourier modes, hence the nonlinear energy transfers from 
these modes may be insignificant. There are large number of Fourier modes associated with 
the lower branch of E (k), and the energy flux possibly results from the nonlinear interactions 
among these modes. For this reason we compare the lower branch of the entropy spectrum 
to either KO or BO scaling. 

Next, we compute the energy fluxes for the velocity and temperature fields for the same 
run. We observe constant fluxes for both the velocity and temperature fields as exhibited in 
Fig. [9j Thus both energy spectra and flux results appear to favor the KO scaling more than 
the BO scaling. Given the kinetic energy spectrum and flux (in the common inertial range), 
we compute Kolmogorov's constant using Eq. (jlj) that yields Kk ~ 2.0 with significant 
error. Considering the uncertainties in the numerical fits, this value is in a reasonable agree- 
ment with Kolmogorov's constant for the fluid or the passive-scalar turbulence measured 
earlier using experiments and numerical simulations. 



We also compute e", e , ka, and fc c using Eqs. ( I7f8|26ll27p . These numbers are listed in 
Table [H The predicted values of e u and e e are in general agreement with the simulation 
results. We observe that k c < k^, which is also evident in the spectra and flux plots. An 
important point to note is that k c ~ 68 is rather large. Hence the arguments presented in 
the earlier subsection for low-P convection will not hold here. More rigorous arguments are 
required to understand the phenomenology for P = 0.2. 

After considering P = 0.2, we turn to convection for P = 1 . 

D. Prandtl number P = 1 

Next we present the energy spectra and fluxes for the kinetic energy and entropy for 
P — 1 at R — 6.6 x 10 6 . Figures fTUl and fTTl exhibit the compensated kinetic energy spectra 
and entropy spectrum respectively. The kinetic energy spectra plots are inconclusive since 
both the compensated plots for the KO and BO scaling are equally flat, albeit at different 
wavenumber ranges. 

The entropy spectrum, shown in Fig. [HJ has two distinct branches similar to low-P cases. 
In agreement with the arguments of Appendix A, the upper branch of the entropy spectrum 
follows E e (2n) ~ n~ 2 . A comparison of the lower curve with the BO or KO scaling indicates 
that neither of the scaling fits well with the numerical data. Figure [12] shows the kinetic 
energy and entropy fluxes along with the compensated kinetic energy flux IP*; 4 / 5 . The flux 
plots are also inconclusive. 

Overall, the numerical results for P = 1 are rather inconclusive. The inverse Bolgiano 
length for P = 1 is approximately 9.0 (see Table 1), hence the phenomenologies predict BO 
scaling for k < 9tt ~ 28 and KO scaling for 28 < k < kd- Clearly the wavenumber range of 
BO or KO scaling is too small to be able to infer any scaling. Also, the arguments put forth 
for the validity of the KO scaling for low-P convection based on the dominance of buoyancy 
force for low wavenumbers cannot be extended to P > 1. 

In the next subsection we discuss the simulation results on convective turbulence for 
P = 6.8. 



E. Prandtl number P = 6.8 



At last, we present the kinetic energy spectrum for P = 6.8 at R — 6.6 x 10 6 . In Fig. [T3] 
we plot the compensated kinetic energy spectra E u (k)k 5/3 (KO) and E u (k)k 11/5 (BO). The 
flat regions in both the plots are rather short, yet the BO line appears to be in a better 
agreement with the numerical results than the KO line. The inverse of Bolgiano length l^ 1 
is around 15.0. Hence according to the convective turbulence phenomenology discussed in 
Section I, the BO scaling should hold for k < vr/^ 1 , and the KO scaling should hold for 
k > ixl'g 1 . The BO scaling appears to be present in our numerical results, but the KO 
scaling is not observable. The dominance of dissipation for modes with k > nl^ 1 in our 512 3 
simulation may be a reason for the absence of the KO scaling. We need higher resolution 
simulation to investigate this issue. 

Figure[T3]exhibits dual branches in the entropy spectrum similar to those discussed earlier. 
The upper spectral curves representing the modes 9(0, 0, 2n) follow k~ 2 scaling as predicted 
in Appendix A. The Fourier modes 0(0,0, 2n) ~ —l/(2nir) as evident from Table HT1 For 
the lower branch, both the KO and BO scaling are not in good agreement with the entropy 
spectrum, yet the BO scaling fits better with the numerical data than the KO scaling. 

Recall that for large-P convection, under the BO scaling, the entropy flux is constant 
but the energy flux varies as /c -4 / 5 (see Eq. (jHJ)). In contrast, in the KO scaling the fluxes 
of the kinetic energy and the entropy are constant. In Fig. [15] we plot both the fluxes as 
well as the compensated kinetic energy flux Il u (k)k 4/5 . We observe that IF(fc) falls rather 
steeply as a function of wavenumbers, but the compensated kinetic energy flux is constant 
in a narrow band of the inertial range. The entropy flux is also a constant for a significantly 
large wavenumber range. Thus the flux results tend to favor BO scaling for P = 6.8. 

Our numerical results on the energy spectra and fluxes are somewhat inconclusive, but 
the BO scaling scaling appears to fit better with the simulation results. 

IV. CONCLUSIONS 

We numerically compute the spectra and fluxes of the velocity and temperature fields 
of convective turbulence using a pseudo-spectral method. We performed these simulations 
for a large range of Prandtl numbers — zero-P, low-P, and large-P. The Rayleigh number 



of our simulation is around a million, which is at the lower end of turbulent convection. 
We apply free-slip and thermal boundary conditions for our simulations. As a result, the 
viscous boundary layer is rather weak, but the thermal boundary layer is quite significant. 
Consequently, our numerical results possibly reflect the scaling for the bulk convective flow. 
The simulation results of kinetic energy and entropy fluxes are in good agreement with their 
estimates using exact relations, thus validating our numerical simulations. 

We find that for nonzero Prandtl numbers, the entropy spectrum shows dual branches. 
Our simulation results indicate that the maximum entropy transfer from the modes 9(n, 0, n) 
and 8(0, n, n) are to the mode #(0, 0, 2n). These observations combined with entropy evolu- 
tion equations yield 9(0, 0, 2n) ~ — l/(2mr) and E e (2n) ~ l/(4n 2 7r 2 ). For large-P convec- 
tion, these predictions fit very well with the upper branch of the entropy spectrum. The 
upper branch however has only a small number of modes, and they probably do not con- 
tribute significantly to the entropy flux. For this reason we compare the lower branch of the 
entropy spectrum to either the Kolmogorov-Obukhov (KO) or the Bolgiano-Obukhov (BO) 
scaling. 

For zero-P convection, the temperature field is active only for very small wavenumbers 
since E e (k) ~ E u (k)/k 4 . Hence, buoyancy is active only for very small wavenumbers leading 
to Kolmogorov's scaling just like in fluid turbulence (E u (k) ~ k~ 5 ^ 3 and Tl u (k) ~ const). 
We observe such behaviour in our numerical simulation. 

The above phenomenological arguments for zero-P convection can be extended to low- 
P convection. For this case, Kolmogorov's diffusive wavenumber k c is much smaller than 
Kolmogorov's dissipation wavenumber kj. Hence, the temperature field will be diffusive 
for k > k c , and the forcing due to buoyancy is active only for low wavenumbers (k < k c ). 
Consequently, we expect Kolmogorov's spectrum for the velocity field for k c < k < k^. 
We numerically compute the energy spectra and fluxes for P = 0.02, and observe diffusive 
spectrum for the temperature field and Kolmogorov's spectrum for the velocity field. Thus 
the phenomenological arguments presented above are in agreement with our simulations. 
For k < k c , the inertial range is too narrow to ascertain any of the KO or the BO scaling. 
A large Rayleigh number simulation could possibly resolve the scaling in this range. 

We have also computed the spectra and fluxes for P = 0.2. For this case, the Kolmogorov- 
Obukhov (KO) scaling appears to fit better than the Bolgiano-Obukhov (BO) scaling with 
the energy spectra and fluxes of the velocity and temperature fields (E u (k) ~ k~ 5 ^ 3 , E 9 ~ 



^-5/ 3 ) Y[ u [k) ~ const, and U 9 (k) ~ const). Numerical results for P = 1 are inconclusive 
regarding the phenomenology. Simulations results for P = 6.8, which is a sample of large-P 
convection, too are inconclusive, however the BO scaling appears to fit better than the KO 
scaling in this case. P = 1 and large-P convection require more refined simulations for 
resolving these issues. 

When we compare our results with earlier experiments and simulations, we observe gen- 
eral agreement with the findings of Cioni et al. 28j| where they reported KO scaling for 
mercury (P = 0.02, low-P). Chilla et al. [2o|, Zhou and Xia |2l|, and Shang and Xia [22] 
performed experiments on water and reported the BO scaling for it. Our simulation re- 
sults are in general agreement with the the above experimental results. A word of caution 
is in order: our simulations use free-slip boundary conditions that differs from the no-slip 
boundary conditions of the experiments. Also, realistic convective flows are quite complex 
due to the presence of boundary layers, anisotropic forcing (buoyancy), plumes, large-scale 
circulation (LSC) etc. all of which have not been analyzed carefully in our analysis. Several 
past numerical simulations and experiments have attempted to study these features 35]. 
Our emphasis in this paper has been on the bulk energy spectrum and fluxes. Note that the 
plumes and LSC typically affect the low-wavenumber regime of the energy spectrum, and 
may not significantly affect the inertial-range isotropic energy spectra being investigated in 
the present paper. 

In summary, we observe the KO scaling for zero-P and low-P convection in our numerical 
simulations. For large-P convection, the numerical results are not very convincing, yet the 
BO scaling matches with the numerical results better than KO scaling. These results are in 
general agreement with some of the earlier experimental and numerical results. We provide 
phenomenological arguments to support KO scaling for low-P and zero-P convection. More 
rigorous theories like renormalization group analysis and very high resolution simulations 
could be very useful in providing further insights into this complex problem. Unfortunately 
convective turbulence simulations beyond 512 3 are prohibitively expensive at this stage. 
Also, more complex features like inhomogeneity, anisotropy, boundary layers need to inves- 
tigated. Future experiments, simulations, and theoretical modeling will hopefully resolve 
this outstanding problem. 



Appendix A: Entropy Spectrum 



The entropy spectrum exhibits dual branches. In this appendix we discuss the reasons 
for this behaviour. We start with the entropy equation for the 6(n,0,n) mode, which is 
d 1^0^ = ^ ^ ^ +^[ U3 (n, 0, n)0*(n, 0, n)\ ^( n 2 n 2 + n 2 k 2 c ) \6(n, 0, n) \ 2 (Al) 

where T e (n, 0, n) is the nonlinear entropy transfer to the mode 6(n, 0,n), and k c = vt/a/2. 
The second term in the RHS is the entropy production rate P e (n,0,n) due to the vertical 
velocity, and the last term provides the dissipation rate of entropy due to thermal diffusivity. 
The entropy equation for the #(0,0, 2n) mode is very similar. We compute T (n, 0,n) and 
P e (n,0,n) from the simulation data, and find these quantities to be highly variable. Yet 
we compute them at a given instant of time in the steady state regime. In this regime, 
d\6(n, 0, n)\ 2 /dt ~ 0, and the dissipation term is also quite small. In Table [TTT1 we list the 
numerical values of T e (n, 0, n), T 9 (0, n, n), P e (n, 0, n), and P e (0, n, n) at an instant. Clearly, 
T e ~ — P e indicating that the entropy generated by is transferred to the higher modes 
by nonlinear transfer. 

From Eq. lAll we can conclude that the 6(n, 0, n) mode gains energy through the entropy 
production term (P e (n, 0, n)), and loses energy to other modes through nonlinear entropy 
transfer (T e (n, 0, n)). When we compute the energy transfers functions explicitly, we find 
that the dominant entropy transfer from 6(n, 0, n) is to the 6(0, 0, 2n). The "mode to mode 



energy transfer" formalism 



40| provides us the entropy transfer rate from 6(0, 0, 2n) to the 



6(n, 0,n) with u(— n, 0,n) acting as a mediator, which is 

5(k|p|q) = -S{2mr(-i)u 3 (-n,0,n)6(n,0,n)6(0,0,2n)} 

= 2mx6(0,0,2n)^[ul(n,0,n)6(n,0,n)} (A2) 

with k = (n, 0,n), p = (0,0, 2n), q = (— n, 0,n). The term stands for the imaginary 
part of the arguments. The above formula has been adopted from the mode-to-mode energy 
transfer formulas for Fourier basis to mixed basis used for the free-slip boundary conditions 
[Eqs. (UMIZD]. 

We compute S'(klplq) using the numerical data at the same instant of time when we 
compute T e (n, 0, n), and compare it with T e (n, 0, n) and the entropy production. As evident 
from the entries of the Table IIHI 

T e (n, 0, n) ~ S(n, 0, n\0, 0, 2n\ - n, 0, n) ~ -P e (n, 0, n). (A3) 



The formulas and relationships for 9(0, n,n) are very similar. The above numerical findings 
indicate that the most dominant entropy transfers to the 9(0, 0, 2n) mode occur from the 
9(n,0,n) and 9(0, n,n). Also, the approximate relationship S(n, 0, n\0, 0, 2n\ — n,0,n) ~ 
— P e (n, 0,n) and the equivalent relationship for the 9(0, n, n) mode yield 

0(O,O,2n)~--i- (A4) 
Irnx 

that matches quite well with the simulation data for P = 6.8 and 0.2 (listed in Table |H| . 
Using the above result we can immediately derive 

E9{Wn) ^4^ (A5) 

that matches well with the upper branch of the entropy spectrum as shown in Figs. [51 El [TTJ 
and [TH The lower branch of the entropy spectrum corresponds to modes other than 
#(0, 0, 2n). The dual branches appear to arise due to the free-slip boundary conditions, and 



they have been observed in the simulation by Vincent and Yuen 



3l| and Paul et al. [32]. 



Note that the dual branches in the entropy spectrum have not been reported for no-slip [34 1 



and periodic boundary conditions 



15 



16]. 

The above arguments that support E e (0,0,2n) ~ 1/n 2 is essentially numerical and phe- 
nomenological that works for P = 6.8 and 1. For lower Prandtl numbers, T e (n,0,n) is not 
approximately equal to S(n, 0, n\0, 0, 2n\ — n, 0, n) possibly due to significant entropy trans- 
fers to other modes, or due to thermal diffusion. Note however that the above quantities are 
within a factor of two, consequently #(0, 0, 2n) ~ holds even for lower Prandtl number 
within a factor of two. 

The dual branches in the entropy spectrum adds complications to the energy fluxes 
discussed in the paper. The temperature modes on the upper branch have significantly 
higher entropy, but they are only a few in numbers. Hence, the nonlinear energy transfers 
arising from the upper branch is possibly insignificant. The number of modes involved in 
the lower branch is quite large, and they are likely to provide the energy flux. 
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FIGURES 




FIG. 1: For R = 6.6 x 10 6 and P = 6.8: (a) the vertical variation of the velocity fluctuations 
averaged over the horizontal planes. The solid, dashed, and dotted lines represent rms values 
of ii£ ms , Uy ms , and u r z ms respectively; (b) the vertical variation of horizontally averaged mean 
temperature. 




FIG. 2: Plot of the compensated kinetic energy spectra E u (k)k 5 / 3 (KO) and E u (k)k u / 5 (BO) vs. 
k for R = 1.97 x 10 4 , P = on 256 3 grid. The DNS spectrum matches with KO spectrum quite 
well. 




FIG. 3: Plot of Kinetic energy flux vs. k for R = 1.97 x 10 4 and P = on 256 3 grid. The constancy 
of the flux in inertial range indicates that zero-P convection follows Kolmogorov's scaling. 




FIG. 4: Plot of the compensated kinetic energy spectra E u (k)k 5 / 3 (KO) and E u (k)k u / 5 (BO) vs. k 
for R = 2.6 x 10 6 , P = 0.02 on 512 3 grid. The KO scaling is in better agreement than BO scaling. 




k 

FIG. 5: Plot of the entropy spectrum for R = 2.6 x 10 6 , P = 0.02 on 512 3 grid. The exponential fit 
in the inset indicates the diffusive nature of the entropy spectrum. The upper part of the entropy 
spectrum corresponds to the 0(0,0, 2n) modes. 




FIG. 6: Plot of the kinetic energy flux (solid line) and the entropy flux (dashed line) vs. k for 
R = 2.6 x 10 6 , P = 0.02 on 512 3 grid. The kinetic energy flux is constant in a narrow inertial 
range indicating agreement with the Kolmogorov's scaling for the velocity field. The entropy flux 
appears to decay rather sharply suggesting diffusive entropy spectrum. 




FIG. 7: Plot of the compensated kinetic energy spectra E u (k)k 5 / 3 (KO) and E u (k)k u / 5 (BO) vs. k 
for R = 6.6 x 10 6 , P = 0.2 on 512 3 grid. The numerical results match better with the KO scaling 
than the BO scaling. 




FIG. 8: Plot of the entropy spectrum vs. k for R = 6.6 x 10 6 , P = 0.2 on 512 3 grid. The lower 
part of the entropy spectrum matches better with the KO scaling than the BO scaling. The upper 
part of the entropy spectrum corresponds to the #(0,0, 2n) modes, and it is in general agreement 
with k~ 2 fit. 




FIG. 9: Plot of the kinetic energy flux (solid line) and the entropy flux (dashed line) vs. k for 
R = 6.6 x 10 6 , P = 0.2 on 512 3 grid. The kinetic energy and entropy fluxes are constant in the 
narrow inertial range indicating a general agreement with the KO scaling. 




FIG. 10: Plot of the compensated kinetic energy spectra E u (k)k 5 / 3 (KO) and E u (k)k u ^ (BO) 
vs. k for R = 6.6 x 10 6 , P = 1 on 512 3 grid. We cannot infer which phenomenology fits better with 
the plots. 




FIG. 11: Plot of the entropy spectrum vs. k for R = 6.6 x 10 6 , P = 1 on 512 3 grid. Both the KO 
and BO lines do not fit with the lower branch of the entropy spectrum. The upper part of the 
entropy spectrum matches with k~ 2 quite well. 




FIG. 12: Plot of the kinetic energy flux (solid line) and the entropy flux (dashed line) vs. k for 
R = 6.6 x 10 6 , P = 1 on 512 3 grid. The dotted line represents U u (k)k A/5 curve. The flux results 
are inconclusive about the nature of scaling. 




FIG. 13: 
k for P 
appears 



Plot of the compensated kinetic energy spectra E u (k)k 5 ^ 3 (KO) and E u (k)k 11 ^ 5 (BO) vs. 
= 6.8, R = 6.6 x 10 6 on 512 3 grid. The fit is somewhat inconclusive, yet the BO scaling 
to fit better with the numerical data than the KO scaling. 




FIG. 14: Plot of entropy spectrum E e (k) vs. k for P = 6.8, R = 6.6 x 10 6 on 512 3 grid. Even 
though both the KO and BO lines do not fit well with the data, BO line is closer to the lower part 
of the spectrum. The upper branch matches with k~ 2 quite well. 




FIG. 15: Plot of the kinetic energy flux (solid line) and the entropy flux (dashed line) vs. k for 
R = 6.6 x 10 6 , P = 6.8 on 512 3 grid. The normalized kinetic energy flux (multiplied by fc 4 / 5 ) is 
also shown in the figure as a dotted line. U u k 4 / 5 and if are constant for a range of wavenumbers. 



TABLES 



TABLE I: Estimates of the viscous dissipation (e u ) and the thermal diffusion rates (e d ) from 
numerical simulation and by using theoretical relationships [Eqs. ([7118])], inverse of the Bolgiano 
length (I]] 1 ), Kolmogorov's dissipative wavenumber (k^), and Kolmogorov's diffusive wavenumber 
(k c ). The reported quantities are nondimensional: e u = (Nu — l)/\/RP, e e = Nu/V RP, kj = 
[(Nu - ljR/P 2 ] 1 / 4 , and k c = [(Nu - ljRP} 1 ^. 



P R Nu e u e u e° e° l^ 1 k d k c 

(estim). (comp.) (estim.) (comp.) 

0.02 2.6 x 10 6 8.5 0.033 0.032 0.037 0.037 5.2 470.3 25.0 

0.2 6.6 x 10 6 17 0.014 0.014 0.015 0.015 8.2 227.0 68.0 

1.0 6.6 x 10 6 32 0.012 0.0082 0.013 0.0085 9.0 108.6 108.6 

6.8 6.6 x 10 6 30 0.004 0.0043 0.0042 0.0042 15.0 44.2 186.2 



TABLE II: Numerical values of 0(0,0,2), 0(0,0,4), 0(0,0,6), and 0(0,0,8) modes for P = 6.8,0.2 
and 0.02. Our phenomenological arguments with numerical ingredients indicate that 0(0, 0,2n) ~ 
-l/(2wr). 



P 0(0,0,2) 0(0,0,4) 0(0,0,6) 0(0,0,8) 

6.8 -0.16 -0.077 -0.050 -0.036 

0.2 -0.15 -0.061 -0.031 -0.017 

0.02 -0.13 -0.040 -0.017 -0.0081 

-0.16 -0.080 -0.053 -0.040 



TABLE III: For high-P (P = 6.8) and low-P (P = 0.2), the numerical values of the nonlinear 
entropy transfer rates T e , entropy production rates P°, and the nonlinear entropy transfer rates 
5(k|p|q) from the mode #(0, 0, 2n) mode to the modes #(n, 0,n) or 6(0, n, n). 





P = 6.8 


P = 0.2 


mode 




5(k|p|q) 


p e 




5(k|p|q) 




(1,0,1) 


-1.2e-4 


-l.le-4 


l.le-4 


-1.3e-3 


-2.7e-3 


2.9e-3 


(0,1,1) 


-1.0e-7 


-7.5e-8 


7.5e-8 


-3.1e-4 


-3.4e-4 


3.6e-4 


(2,0,2) 


-7.0e-7 


-5.0e-7 


5.5e-7 


-2.5e-6 


-2.6e-5 


3.4e-5 


(0,2,2) 


-1.6e-6 


-1.2e-6 


1.3e-6 


-6.1e-5 


-4.9e-5 


6.3e-5 


(3,0,3) 


3.0e-7 


2.1e-7 


-2.2e-7 


-7.0e-7 


-1.4e-6 


2.3e-6 


(0,3,3) 


-1.0e-7 


-9.3e-8 


1.0e-7 


-5.7e-6 


-5.3e-6 


9.2e-6 


(4,0,4) 


3.0e-7 


-3.2e-7 


3.5e-7 


-8.0e-7 


-1.0e-6 


2.4e-6 


(0,4,4) 


-4.0e-7 


-4.0e-7 


3.9e-7 


1.0e-6 


-7.0e-7 


-1.6e-6 



